Project Euler tutorial

For those who work on the project Euler problems, this article contains a spoiler for bonus problem n°18i. Though, you are probably not looking for the kind of solution explained here

The goal of this article is to provide a tutorial on how to solve problems from project Euler without doing a lot of maths, by using Télécom's computers. I'm going to take this problem as an example, and explain my answer:

While there is clever reasonning to simplify the problem, the most straightforward way is the brute-force manner.

First, we compute $R(p)$: It is possible to do $p/2$ iterations instead of $p$, by gathering the $x$ and $-x$ members of the product:

By the way: we don't try to optimize in a way "$O(n\log n)$ instead of $O(n^2)$ through a clever algorithm", but we still try to gain a factor of $2$ whenever we can, as long as it requires little to no brain-work.

We can rewrite the previous result to: $-((x·(x^2-3){)}^2-16)$. For $x=0$, the polynomial evaluates to $4$. As $p$ is odd, there is no issue such as $p/2\equiv -p/2$ Thus, $R(p)=4\times {\prod }_{x=1}^{\lfloor p/2\rfloor }-(x·(x^2-3)-16)$ As values are quite large (over ${10}^9$), we need to apply a modulo after each multiplication to avoid overflows. Finally, whenever we reach a value of $0$, we directly return without computing the other polynomials. Thus, we obtain the following (Rust 🦀) code:

fn R(p: i64) -> i64 {
    let mut prod = 4; // polynomial for x = 0 gives 4
    for x in 1..=p / 2 {
        let poly = x * ((x * x - 3) % p) % p;
        let poly = (poly * poly - 16) % p;
        prod = (- prod * poly) % p;
        if prod == 0 {
            return prod;
        }
    }
    (p + prod) % p // So that the result stays non-negative
}

Then, we compute $R(p)$ for each prime $1000000000\<p\<1100000000$ using multithreading to speed up the process, and do the sum:

const INF: usize = 1_000_000_000;
const SUP: usize = 1_100_000_000;
fn bonus() {
    assert_eq!(R(11), 0);
    assert_eq!(R(29), 13);
    let primes_list = Arc::new(primes_list_until(SUP)); // The primes_list_until function was home-made
    let nb_threads = 64;
    let mut handle = Vec::new();
    for thread_id in 0..nb_threads {
        let primes_list = primes_list.clone();
        handle.push(thread::spawn(move || {
            let mut res = 0;
            for &p in primes_list.iter().skip(thread_id).step_by(nb_threads) {
                if p > INF {
                    res += R(p as i64);
                }
            }
            res
        }))
    }
    let mut res = 0;
    for t in handle.into_iter() {
        res += t.join().unwrap();
    }
    println!("{res}");
    return;
}

The main point is this line: let nb_threads = 64;. By looking at this webpage: https://lames.telecom-paris.fr/ (you may need to be connected to the Campus-Telecom WiFi, or using the VPN), you can find some machines such as lame25 or lamedell18 which are available for students, and have 128 threads. Though, take care not to monopolize a machine: they are used by teachers and PhD students at Télécom, so use at most a third or half of their CPU power.

Finally, after waiting for a bit (4 to 5 days), you get the answer. And while you used much more electricity and time than the user pacome_f, who spent "18 seconds with pypy", you probably saved a lot of brainpower compared to him!